Excitations in disordered bosonic optical lattices 
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Spectral excitations of ultracold gases of bosonic atoms trapped in one-dimensional optical lattices 
with disorder are investigated by means of the variational cluster approach applied to the Bose- 
Hubbard model. Qualitatively different disorder distributions typically employed in experiments 
are considered. The computed spectra exhibit a strong dependence on the shape of the disorder 
distribution and the disorder strength. We compare alternative results for the Mott gap obtained 
from its formal definition and from the minimum peak distance, which is the quantity available from 
experiments. 
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Interacting many-body systems with disorder are fas- 
cinating and challenging from both the experimental as 
well as the theoretical point of view. Understanding 
disordered bosonic systems has been of great interest 
ever since the pioneering works on the Bose-Hubbard 
(BH) model [1|, which describes strongly interacting lat- 
tice bosons. Originally, the disordered BH model has 
been used to approximately describe various condensed 
matter systems, such as superfluid helium absorbed in 
porous media @, Q, superfluid films on substrates [J], 
and Josephson junction arrays [5[. However, seminal ex- 
periments on ultracold gases of atoms trapped in optical 
lattices shed new light on interacting bosonic many-body 
systems, as these systems provide a direct experimental 
realization of the BH model 0, 0|- Intriguingly, these 
experiments allow to observe quantum many-body phe- 
nomena, such as the quantum phase transition from a 
superfluid to a Mott state Q ■ The condensate of atoms 
can be driven across this phase transition by gradually 
increasing the laser beam intensity, which is directly re- 
lated to the depth of the potential wells. There is a large 
experimental control over the system parameters such as 
the particle number or lattice depth, and in addition the 
parameters are tuneable over a wide range. While opti- 
cal lattices provide a very clean experimental realization 
of strongly correlated lattice bosons, they can be used to 
study disordered systems on a very high level of control as 
well. Disorder can be added to the regular optical lattice 
by several techniques, such as by superposing additional 
optical lattices with shifted wavelength and beam angles 
[2-Q3, laser speckle fields |15l - [l8j . or including atoms of 
a different species acting as impurities [19| . 

The disordered BH model has been widely investigated 
in the literature. Most of the work has been devoted to 
study the phase transitions occurring in the disordered 
BH model [20l - l35| . At zero temperature the ground-state 
phase diagram depending on the chemical potential, the 
tunneling probability of the particles and the disorder 
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FIG. 1. (Color online) The periodic optical lattice potential is 
locally modified by disorder. In the illustration the disorder 
is bounded by e* corresponding to a situation obtained by the 
superposition of two incommensurate optical lattices. 



strength, consists of three different phases: the Mott in- 
sulating, the superfluid, and the Bose glass phase. The 
first two phases are already present in the pure BH model, 
while the latter is a distinct feature of the disordered sys- 
tem. The Bose glass phase is characterized by being gap- 
less and compressible, however, due to disorder the phase 
coherence does not extend over the entire system in con- 
trast to the superfluid phase. While the phase diagram of 
the disordered BH model has been extensively studied, 
up to now spectral properties have not yet been inves- 
tigated theoretically, even though they are experimen- 
tally accessible by Bragg spectroscopy, which allows to 
extract the wave vector dependent excitation energies of 
the system [3a - l38| . For low-energy excitations (typically 
uj/2-k < 10kHz) the response of the Bragg spectroscopy 
is ascribed to the structure factor, which corresponds to 
density fluctuations, whereas, for high excitation ener- 
gies (u)/2ir > 30kHz) the single-particle spectral function 
can be extracted [39|, Eoj]. In this case, atoms excited 
into high-energy bands are not expected to interact with 
the lower bands. Therefore, the measured spectrum de- 
scribes a convolution (without vertex corrections) of the 
density of states of the lower occupied with the upper un- 
occupied band. From this convolution the spectral func- 
tion of the lower strong-correlated bands can be extracted 

In the present paper we study in detail the spectral 
properties of bosonic atoms in a disordered optical lat- 
tice, modeled by the one-dimensional BH model in which, 
for simplicity, the harmonic trap potential has been ne- 
glected. We focus on the strongly correlated regime and 



evaluate spectral properties for the hopping strength to 
the on-site interaction ratio of size t/U = 0.05. This ratio 
corresponds to an optical lattice depth of 11 En , where 
Er is the recoil energy [41(. To evaluate this depth we 
considered laser beams with wave length A = 830 nm and 
the scattering length of rubidium a s (=a 5 nm. In particu- 
lar, we focus on the similarities and differences between 
the two common experimental methods used to induce 
disorder: (i) superposition of two laser fields with in- 
commensurate wavevectors on the one hand, and (ii) the 
addition of a laser speckle field on the other hand. The 
disordered BH model and the disorder distributions gen- 
erated by the two experimental methods are introduced 
in Sec. HI 

We investigate the spectral properties of the disordered 
BH model numerically using the variational cluster ap- 
proach [42| ■ In this work we present the first application 
of the variational cluster approach on disordered, inter- 
acting many-body systems. A description of the method 
is given in Sec. |TIJ Our results on the excitations of the 
disordered BH model are presented in Sec. IIHI The cru- 
cial point about our investigations is that the two experi- 
mental approaches described previously produce disorder 
with rather different distributions leading, in turn, to dif- 
ferent distributions of the spectral weight, as we show in 
Figs. ([2]) and ([5]). In addition, we present two results for 
the Mott gap. The first one, termed A, is obtained from 
the smallest excitation energy. This is, in principle, the 
formally correct version, but it is hard to determine ex- 
perimentally. The second one, termed A mp , is evaluated 
from the minimum peak distance in the spectral weight, 
which is the accessible quantity in experiments. The two 
quantities exhibit a different behavior as shown in Fig. [4] 
Finally, in Sec. lIVl we conclude and summarize our work. 



is distributed according to a given probability distribu- 
tion function (pdf) p(rj), which has to satisfy the nor- 
malization condition J p{n) dn = 1. The average of a 
quantity A,, with respect to the pdf p(j]) is given by 
X p = (X) p = J p{rj) X v dn. We consider e 2 ; as iden- 
tically and independently distributed random variables 
leading to p(rf) = Y\ i=1 q{u), where q(ei) is a pdf de- 
scribing the disorder generated in the experiment. Ex- 
perimentally, there are two main ways of introducing dis- 
order, corresponding to two different pdf's. On the one 
hand, superposing two incommensurate optical lattices 
yields a shifted j3 distribution as demonstrated in Ap- 
pendix [5J The distribution is bounded by a maximum 
disorder strength e*. This approach is used, for exam- 
ple, in the setup by Fallani et al. [lflf , where two laser 
beams with wave lengths Ai = 830 nm and A2 = 1 076 
nm are superposed. Since these two wave lengths are 
highly incommensurable the disorder can be regarded as 
truly random, see Ref. y. Alternatively, disorder can be 
generated by superposing a laser speckle field [4j, [44j to 
the regular optical lattice. In this case, one obtains an 
exponential distribution q{ti) = flfe) exp(— ej/e)/e where 
e specifies the mean disorder strength |43|,|44| . In our cal- 
culations the distribution is additionally shifted about its 
median e In 2 to avoid modifications of the chemical po- 
tential due to disorder. The exponential distribution is 
asymmetric and unbounded and thus, strictly speaking, 
the Mott phase, which is controlled by the extrema of 
the distribution [23|, [2J, [3l], |32[ does not occur anymore. 
However, the goal of the present paper is to mimic the 
results of the experimental measurements, for which the 
effects of the tail of the distribution are too small to be 
observable. Therefore we introduce a cut off for the dis- 
order distribution of 4 1 1n 2. In other words, realizations 
far off the median are not considered. 



I. DISORDERED BOSE-HUBBARD MODEL 

The Hamiltonian of the BH model with on-site disorder 
[l| is given by 

H= -tJ2 blbj + ^Y.^ fa - l)+J2&-^ > (^ 
(*, J) l * 

where the operators b\ (6^) create (annihilate) bosonic 
particles at lattice site i. The first sum in Eq. ((T|) is re- 
stricted to nearest-neighbor sites. The hopping strength 
is denoted as t, U is the on-site repulsion and u stands 
for the chemical potential, which controls the total par- 
ticle number N p = J2i n i = J2i K &, ■ The site-dependent 
random variable e» introduces disorder. In the following 
calculations and discussions we use the local interaction 
U as unit of energy. 

One sample-configuration of disorder is denoted as 
77 = (ei, €2 ■ ■ ■ ejv), where N is the number of lattice 
sites. The disorder modifies locally the depth of the op- 
tical lattice, see Fig. Q] for illustration. Random disorder 



II. VARIATIONAL CLUSTER APPROACH 

We employ the variational cluster approach (VCA) to 
investigate the phase boundary and spectral properties 
of the disordered BH model. VCA has been originally 
proposed for fermionic systems by M. Potthoff et al. in 
Ref. I42t and has been extended to bosonic systems in 
Refs. [45H47| . The mathematical concept for the treat- 
ment of disordered systems by means of VCA can be 
found in Ref. [4g. In this paper we present the first ap- 
plication of VCA for disordered systems and on top of 
that, as a new contribution to VCA theory, we extend 
the so-called Q-matrix formalism [47|, |49|, [5(| to a formu- 
lation suitable for disordered systems. 

The basic concept of VCA is that the exact self-energy 
H ex of a certain interacting system is given by the sad- 
dle point of a functional fl[£l, which becomes the exact 
grand potential at Y. = H ex [5ll.l52|. The important point 
is that the functional Q can be computed exactly for a 
restricted set of self-energies associated with a so-called 
reference system, described by an Hamiltonian H' , which 
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FIG. 2. (Color online) Spectral function A p (k, u) for disorder generated by superposition of incommensurate optical lattices. 
The parameters are t = 0.05, fi = 0.45, and (a) e* = 0.0, (b) e* = 0.1, (c) e* = 0.2, and (d) e* = 0.3. 



shares its interaction part with the physical system. Typ- 
ically, the reference system H' is a decomposition of the 
physical system into independent clusters of a given size 
L, so that its Green's function can be determined numer- 
ically. 

Disorder is treated by introducing a reference system 
which shares, in addition to the interaction part, the dis- 
order distribution with the physical system. The expres- 
sion for the averaged grand-potential within VCA then 
reads 0: 



- ' / 



fii + Trbif-G-^-Trml 



-g;- 1 ] 



(2) 



where £l' p and G' p are the exact disorder averaged grand 
potential and Green's function of the reference system, 
respectively, which can be easily evaluated numerically 
with high accuracy. The averaged Green's function of 



the physical system reads G p 



(g;- 1 



v) 



where 



V = G ~ — Gq , and Go and G are the noninteract- 
ing Green's functions of the pure physical and reference 
systems, respectively. 

In Ref. [47| it has been discussed for a pure system that 
it is expedient for the evaluation of the traces in Eq. ([2]) 
to use the Q-matrix formalism. This formalism is based 



on the Lehmann representation of the Green's function 
G' of the reference system and that of the physical sys- 
tem G. In the present paper, we extend the Q-matrix 
formalism to the case of disorder. For a specific disorder 
configuration rj we have 



G; = Q„ g ;sQ+, 



(3) 



where g' = u — A v and (/\. v ) rr ' = AJ? S rr i are the poles 
of the reference Green's function, see Ref. [47( for details. 
In praxis a specific number M of disorder configurations 
77 (each consisting of L disorder realizations ti) has to be 
sampled to compute the average of a quantity X v leading 
to X p = (X) p w l/ME Q I, a . The averaged Green's 
function G' thus reads 
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To exploit the Q-matrix formalism we write the averaged 
Green's function in a form similar to Eq. ([3]) and define 
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FIG. 3. (Color online) (a) Density of states N p (lj) evaluated 
for the same parameters as in Fig. [2] (b) Illustration of the two 
"gaps" A and A mp within a blowup of the the data for disor- 
der strength e* = 0.3. According to our definition (see text), 
A corresponds to the formal definition of the Mott gap based 
on minimal excitation energies, whereas A mp corresponds to 
the minimum peak distance in the spectral weight, which is 
the quantity available from experiments. 



where Q = (Q V1 /VM, Q, l2 /VM . . .), S = diag(S, S, . . .) 
and g' = diag(g' , g' V2 , . . .). With that we can proceed in 
the same way as for pure systems and write the averaged 
Green's function of the physical system as 



G P = Q 



1 



-(A-SQtVQ) 



SQt 



where A = diag(A m , l\ m . . .). By diagonalizing the ma- 
trix M = A - SQ f VQ = XDX- 1 wc obtain the poles 
D of the physical Green's function and are thus able to 
evaluate the grand potential fi p (x). 




FIG. 4. (Color online) Gap A obtained from (a) the smallest 
excitation energies and (b) A mp obtained from the minimum 
peak energy difference of the occupied and unoccupied bands. 



III. RESULTS 

We first investigate spectral properties for disorder re- 
alizations sampled from the shifted /? distribution, which 
is realized experimentally by superposing two incommen- 
surate optical lattices, see Appendix lAl for details on the 
distribution. 

In particular, spectral functions A p (k, w) are evalu- 
ated for hopping strength t = 0.05, chemical poten- 
tial ii = 0.45 and various strengths of disorder e* = 
{0.0, 0.1, 0.2, 0.3}, see Fig.H The corresponding den- 
sities of states N p (uj) are shown in Fig. [3] 

For the numerical evaluation by means of VCA we use 
AI = 256 disorder configurations, a reference system of 
size L = 8, and the VCA parameters x = {/i, t, S}. The 
parameter 5 is an additional on-site energy located at the 
boundaries of the cluster, whose introduction drastically 
improves the results, as shown in Ref. [53J , The artifi- 
cial broadening parameter is chosen to be + = 0.05 for 
spectral functions and + = 0.01 for densities of states. 



As mentioned in the introduction the occupied part of 
the single-particle spectral function is experimentally ac- 
cessible by high-energy Bragg spectroscopy. For the pure 
system we obtain the well-known cosinelike shaped bands 
reminiscent of the dispersion of free particles on a lattice. 
For increasing disorder strength e* the minimal gap be- 
tween the occupied and the unoccupied band shrinks. 
In addition, the bands become broader as the poles of 
the Green's function for a specific wave vector k are dis- 
tributed over a large energy range. This behavior can 
also be seen in the density of states. Furthermore, for 
large disorder the bands seem to split in two sections 
separated by a pseudogap around k = it. This is a pe- 
culiarity of this disorder configuration which to some ex- 
tend resembles a binary distribution. An ordered binary 
distribution would double the unit cell thus producing a 
true gap at k = n. 

Experimentally the gap A present in Mott phase can 
be determined for instance by lattice modulation [54| or 
by Bragg spectroscopy [361 ]. In these experiments the 
amount of energy transferred to the system is related to 
the width of the central peak observed in time of flight 
images [54J . The width of the central peak is measured 
for various energies leading to the excitation spectrum. 
Indeed for increasing disorder strength e* a broadening 
of the excitation band has been observed [10| which is 
qualitatively in agreement with our results for the spec- 
tral function. In addition, it could be shown experimen- 
tally that the weight moves to lower excitation energies 
for increasing e* [10] - However, it is rather difficult to 
extract the precise value of the gap from the measured 
excitation spectra. Strictly speaking, the gap A is de- 
fined as the energy difference of the lowest lying poles in 
the occupied and unoccupied bands, respectively. Yet, 
our results show, that these poles quite generally carry 
only very little spectral weight. This means that it is vir- 
tually impossible to detect them in the experiment. For 
this reason, it might be useful to introduce the notion of 
a gap A mp , that corresponds to the experimental situa- 
tion and is determined by the distance between the max- 
ima of the spectral weight observed in the center of the 
Brillouin zone. In Fig. @] we compare A [Fig.@](a)] with 
A rap [Fig. 2](b)] which is obviously always larger than A. 
For increasing hopping strength t the gap A mp decreases, 
however, it still remains finite (i. e., A mp sa 0.2) at values 
of t for which the gap A determined from the smallest 
excitation energies is already almost zero. Therefore, a 
peak at finite energy will be observed in the experimental 
data, even when the system is already in the Bose-glass 
phase. Additionally, it is important to mention that the 
smallest excitation energy gap A can be predicted ana- 
lytically for systems with infinitely many disorder realiza- 
tions. In this case A is controlled by the maximum disor- 
der strength e* only. In particular, the phase boundary 
is shifted by ±e* as there exist always rare regions where 
the chemical potential is either increased or decreased by 
e* [23, [24j. Here, however, we take into account a finite 
number of L M = 2 048 random values for the on-site 
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FIG. 5. (Color online) Spectral properties for a disorder dis- 
tribution generated by laser speckle fields. The parameters 
are t — 0.05, and fi = 0.35. Panel (a) shows the spec- 
tral function A p (k, u) evaluated for e = 0.15 and panel (b) 
shows the densities of states for various disorder strengths 
e= {0.0, 0.05, 0.1, 0.15}. 



energies e,; leading to larger gaps, since it is very unlikely 
that all values a of one realization r\ are close to the 
extreme cases ±e*. This resembles more closely the ex- 
perimental situation in which only a limited number of 
disorder realizations can be detected. 

The second kind of disorder we are addressing in this 
paper follows the shifted exponential distribution, which 
is generated by superposing a laser speckle field. Spectral 
properties of the Bose-Hubbard model for this disorder 
distribution are shown in Fig.[SJ In contrast to the disor- 
der realized by the superposition of two incommensurate 
optical lattices, the spectral signatures evaluated for the 
exponential distribution clearly exhibits an asymmetric 
shape. This becomes particular visible when comparing 
the densities of states for various disorder strengths, see 



Fig.[5](b). In particular, the poles are smeared out to- 
ward higher excitation energies and thus reflect the tail 
of the exponential distribution. No pseudogap behav- 
ior around k = n is observed here, due to the shape of 
the exponential distribution which is in contrast to the j3 
distribution not peaked at the edges. 
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IV. CONCLUSIONS 
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In the present work we investigated for the first 
time spectral properties of the disordered Bose-Hubbard 
model. In view of a realistic description of the experimen- 
tal results, we focused on disorder distributions which 
are relevant for ultracold gases of atoms in optical lat- 
tices. In particular, we studied the differences between 
disorder realized by the superposition of incommensu- 
rate laser fields with the one obtained by laser speckle 
fields. In both cases we evaluated spectral functions and 
densities of states and showed that the resulting spec- 
tral weight strongly depends on the underlying shape of 
the disorder distribution. Furthermore, we determined 
the gap present in the Mott phase for disorder generated 
by incommensurate optical lattices. On the one hand, we 
evaluated the gap A from the minimal excitation energies 
of the system and on the other hand we determined the 
gap A mp from the minimum peak distance in the spec- 
tral weight located at the center of the Brillouin zone. 
Whereas A cannot be observed directly in the experiment 
since the low-energy excitations carry very little spectral 
weight, A mp is directly measurable. Furthermore, A mp 
is always larger than A and thus A mp remains finite even 
at the Mott to Bose-glass transition. In our calculations, 
we neglected the harmonic trap potential present in the 
experiments with ultracold gases of atoms. In principle, 
this effect can be included in our formalism, however, 
with a significantly major effort which goes beyond the 
goal of the present work. There are two cases in which 
neglecting the trap potential is justified. As in the or- 
dered case, one can expect for sufficiently small disorder 
strength the trap potential to lead to multiple ringlikc 
regions which are alternately Mott gapped and gapless. 
Quite generally, one should be able to choose the param- 
eters so that the volume of the gapless regions is much 
smaller than the one of the N p = N Mott region in which 
we are interested. In this case, the spectrum will display 
a nonvanishing weight within the gap originating from 
the gapless regions. However, we expect this to be small 
enough, at least for lj ^ 0, so that the peaks defining 
the experimental gap A mp remain discernible. Alterna- 
tively, spectroscopy experiments probing the system lo- 
cally should be able to probe directly the Mott insulating 
region with no contributions from the gapless areas. 



FIG. 6. (Color online) Total lattice potential created by the 
superposition of the main optical lattice with depth si and 
the disorder lattice with depth S2- 
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Appendix A: Disorder distribution generated by 
incommensurate optical lattices 

In this Appendix we show that the potential distribu- 
tion, which results from the superposition of two optical 
lattices, follows a shifted (3 distribution. In particular, 
we focus on the experimental setup of Fallani et al. [lCf , 
who used for the main optical lattice a laser at wave- 
length Ai = 830 ran. Disorder is generated by superpos- 
ing an additional lattice created from a weak laser beam 
at A2 = 1 076 nm. The resulting potential is given by 

V(x) = si -Eri sin 2ttx/\i + S2 -Er2 sin 2ttx/\2 , 

where x is the spatial position, si and S2 are related 
to the depth of the potential generated from the laser 
beams at Ai and A2, respectively The constants Em 
and Er2 are the corresponding recoil energies. In this 
Appendix, the lattice depths si and S2 will always be 
stated in units of their recoil energies. The depth S2 of 
the disorder-inducing wave is related to the maximal dis- 
order strength e* by e* = S2/2. Since the wave lengths Ai 
and A2 are incommensurable the disorder imitates a true 
random behavior [9(. Here we reconstruct the histogram 
of the disorder distribution and find a mathematical ex- 
pression, which reproduces the behavior of these physical 
systems. 

In the experiments typical values for the lattice depths 
are si = 16 and S2 = 2, see Ref. [I0|. Figure |5] compares 
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FIG. 7. (Color online) Distribution of the on-site energy e; ob- 
served in the experiment (histogram) compared with random 
samples drawn from the shifted j3 distribution (solid line). 



the potential V{x) for the previously mentioned values of 
,si = 16 and si =1 with the pure case, where the second 
laser beam at A2 is switched off, i. e., s\ = 16 and s<i = 0. 
It can be seen that the on-site energy varies for distinct 
lattice sites mimicking a random potential. Actually, the 
parameter S2 just scales the disorder and thus our con- 
siderations are valid for arbitrary disorder strength e*. 
To evaluate a histogram of the energy distribution we 
subtract the disordered potential at the lattice sites from 
the pure potential and shift the difference by its mean 
e*. This yields a distribution which is centered around 
zero. The shift could have been absorbed as well in the 
definition of the chemical potential fi. The histogram for 
2 048 lattice sites is shown in Fig. [7] In the center the 
distribution is rather flat, yet, there are important fea- 



tures at the boundary. Such a distribution can be well 
described by a shifted ft distribution. 

The j3 distribution qp(u\a, b), which is defined on the 
interval u € [0, 1], is given by 



qp(u\a, b) 



B(a, b) 



u a - x (1 - u) 



6-1 



where 



B(a, b)= f dpp a - 1 (l-pf- 1 

Jo 

The expectation value of qp{v\a, b) is 

, . a 



r(#) 

T(a + b) 



a + b 



and its variance is 



var(u) 



W(l-(")) 
a + b + 1 



(Al) 



(A2) 



To obtain a distribution which is symmetric around zero 
and bounded by [— 1, 1] we set (u) = 0.5, shift the whole 
P distribution by this value and scale it by a factor of 
2. Resolving Eqs. (|A1[) and (|A2[) under the condition 
that (u) = 0.5 leads to a = b = [var(u)/4 - l]/2. For 
var(u) > 1/12 (where 1/12 is the variance of the uni- 
form distribution), the probability density is shifted to- 
ward the boundaries of the distribution. In particular we 
set var(u) = 0.12. Finally, we draw 2 048 samples from 
this specific shifted /3 distribution, which we denote as 
q(ei/e*). The resulting distribution of the 2 048 samples 
is indicated by the solid line in Fig. which reproduces 
well the distribution obtained from the experiment with 
two incommensurate optical lattices. 
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